#AUTHOR:YEJunping
#此程序是为了以预处理之后的All_sample_merged.HC.VQSR.genotype.heter(another_all.txt)为标准，比对$my_dir目录下的各个文件
#T为在标准文件中能匹配上的，F为匹配不上的，N为没有找到该位点(如果T,F都为0，N有数值，表示该样本不存在标准文件中)
#!/usr/bin/perl
$all="/mnt/md1200/6/yjp/ASM/another_all.txt";
open(GENE,"$all");
%M18=();
%M20=();
%M22=();
%M40=();
%M26=();
%M27=();
%M5=();
%M36=();
%M30=();
%M47=();
%M49=();
%M7=();
%X1T=();
$line=<GENE>;
$line=<GENE>;
$M18{'key'}="TRUE";$M20{'key'}="TRUE";$M22{'key'}="TRUE";$M40{'key'}="TRUE";$M26{'key'}="TRUE";$M27{'key'}="TRUE";$M5{'key'}="TRUE";
$M36{'key'}="TRUE";$M30{'key'}="TRUE";$M47{'key'}="TRUE";$M49{'key'}="TRUE";
$M7{'key'}="TRUE";$X1T{'key'}="TRUE";
while ($line ne ""){
chomp($line);
@arr=split(/\t/,$line);
if(length($arr[4])==3){$M18{$arr[1]}=$arr[4];}
if(length($arr[5])==3){$M20{$arr[1]}=$arr[5];}
if(length($arr[6])==3){$M22{$arr[1]}=$arr[6];}
if(length($arr[7])==3){$M40{$arr[1]}=$arr[7];}
if(length($arr[8])==3){$M26{$arr[1]}=$arr[8];}
if(length($arr[9])==3){$M27{$arr[1]}=$arr[9];}
if(length($arr[10])==3){$M5{$arr[1]}=$arr[10];}
if(length($arr[11])==3){$M36{$arr[1]}=$arr[11];}
if(length($arr[13])==3){$M30{$arr[1]}=$arr[13];}
if(length($arr[15])==3){$M47{$arr[1]}=$arr[15];}
if(length($arr[16])==3){$M49{$arr[1]}=$arr[16];}
if(length($arr[18])==3){$M7{$arr[1]}=$arr[18];}
if(length($arr[19])==3){$X1T{$arr[1]}=$arr[19];}
$line=<GENE>;
}
close(GENE);

$my_dir='/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/ASM/varscan';
opendir(DIR, $my_dir);		#打开目录
@dirs = grep {/^*.snp/} readdir(DIR);
open (OUTRESULT,">>result.log.txt");
print OUTRESULT "filename\tT\tF\tN\n";
close(OUTRESULT);
foreach $d (@dirs){
if($d eq '.' || $d eq '..'){
next;
}
print "$d is operating\n";
open (OUTRESULT,">>result.log.txt");
$newdir = $my_dir."/$d";
@id_pre=split(/\./,$d);
@id=split(/_/,$id_pre[0]);	#文件名,如果要修改程序以适应其他用法，这里修改文件名的提取很重要！

if (exists ${$id[0]}{key}){
$N=$T=$F=0;
open(SAMPLE,"$newdir");
open (OUTSAMPLE,">$d");
$sample=<SAMPLE>;
chomp($sample);
$header=$sample."\texits\n";
print OUTSAMPLE "$header";

$sample=<SAMPLE>;
while($sample ne ''){
	chomp($sample);
	@arrary=split(/\t/,$sample);
	$sample_id=join ":",$arrary[0],$arrary[1];	#需要比对的唯一id
	$sample_genotype=$arrary[2].$arrary[3];		#需要比对的基因型
	if (exists ${$id[0]}{$sample_id}){			#哈希表里是否存在这个id
	if (${$id[0]}{$sample_id}=~/$sample_genotype/){$merge_outline=$sample."\tT\n";$T++;print OUTSAMPLE "$merge_outline";}
	else {$merge_outline=$sample."\tF\n";$F++;print OUTSAMPLE "$merge_outline";}
	}
	else {$merge_outline=$sample."\tN\n";$N++;print OUTSAMPLE "$merge_outline";}
	#print"$sample_id\t$sample_genotype\n";
	$sample=<SAMPLE>;
}
print OUTRESULT "$d:\t$T\t$F\t$N\n";
close(SAMPLE);
close(OUTSAMPLE);
close(OUTRESULT);
}
else {
$N=$T=$F=0;
open(SAMPLE,"$newdir");
open (OUTSAMPLE,">$d");
$sample=<SAMPLE>;
chomp($sample);
$header=$sample."\texits\n";
print OUTSAMPLE "$header";
$sample=<SAMPLE>;
while($sample ne ''){
chomp($sample);
	@arrary=split(/\t/,$sample);
	$sample_id=join ":",$arrary[0],$arrary[1];	#需要比对的唯一id
	$sample_genotype=$arrary[2].$arrary[3];		#需要比对的基因型
	if (exists ${$id[1]}{$sample_id}){			#哈希表里是否存在这个id
	if (${$id[1]}{$sample_id}=~/$sample_genotype/){$merge_outline=$sample."\tT\n";$T++;print OUTSAMPLE "$merge_outline";}
	else {$merge_outline=$sample."\tF\n";$F++;print OUTSAMPLE "$merge_outline";}
	}
	else {$merge_outline=$sample."\tN\n";$N++;print OUTSAMPLE "$merge_outline";}
	#print"$sample_id\t$sample_genotype\n";
	$sample=<SAMPLE>;
}
print OUTRESULT "$d:\t$T\t$F\t$N\n";
close(SAMPLE);
close(OUTSAMPLE);
close(OUTRESULT);
}

}